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FORECASTING  WATER  TEMPERATURE 
DECLINE  AND  FREEZE-UP  IN  RIVERS 


H.T.  Shen,  E.P.  Foltyn  and  S.F.  Daly 


INTRODUCTION 

The  formation  of  an  ice  cover,  or  freeze-up,  in  a  river  is  an  important  phenomenon  that  affects 
winter  operations.  The  ability  to  provide  accurate,  long-range  forecasts  in  the  fall  for  the  date  of 
initial  ice  formation  in  a  river  is  important  in  scheduling  the  shipping  season,  operating  hydro- 
lectric  plants,  planning  flow  regulation,  and  scheduling  the  installation  of  ice  control  structures 
and  devices.  River  ice  forms  when  the  water  temperature  declines  from  its  summer  maximum  to 
freezing.  The  decline  of  the  water  temperature  is  governed  by  the  heat  exchange  at  the  air/water 
interface. 

By  assuming  that  the  heat  exchange  is  directly  proportional  to  the  temperature  difference  be¬ 
tween  the  water  surface  and  the  air,  Rodhe  (1952)  developed  an  iterative  relationship  between 
the  mean  daily  air  temperature  and  ice  formation  in  the  Baltic  Sea.  Using  Rodhe’s  formulation, 
BileOo  (1964)  developed  a  method  for  predicting  river  and  lake  ice  formation.  Greene  (1983) 
recently  developed  a  procedure  for  forecasting  freeze-up  in  the  St.  Marys  River  by  introducing  air 
temperature  forecasts  into  BUello’s  method.  The  Rodhe-BileDo  method  considers  the  decay  of 
water  temperature  as  a  function  of  a  single  site-specific  constant.  Their  formulation  does  not  in¬ 
clude  die  convective  effects  of  the  river  flow  and  is  therefore  not  strictly  valid  for  rivers.  Poulin 
et  al.  (1971 )  developed  a  probability  forecast  method  for  the  water  temperature  and  freeze-up 
dates  in  the  St.  Lawrence  River.  Instead  of  incorporating  air  temperature  forecasts  in  their 
method,  Poulin  et  al.  determined  air  temperature  regimes  in  terms  of  probabilities  from  the  past 


record,  and  they  forecast  freeze-up  for  different  probable  air  temperature  regimes.  They  used  the 
normal  decline  for  die  water  temperature  at  the  upstream  end  of  the  river.  Adams  (1976)  and 
Assel  (1976, 1977)  developed  freeze-up  forecast  models  for  die  upper  St.  Lawrence  River  based 
on  analytical  formulations  similar  to  that  of  Poulin  et  al.  (1971 ).  Assel  incorporated  air  tempera¬ 
ture  forecasts  in  this  model,  but  die  normal  air  temperature  was  assumed  to  be  constant  during 
each  half-month  period. 

In  this  report  the  normal  pattern  of  air  temperature  variation  during  a  year  is  examined,  and  the 
relationship  between  variations  in  air  temperature  and  water  temperature  is  analyzed.  On  the  basis 
of  these  analyses,  an  improved  analytical  method  for  predicting  water  temperature  decline  and 
freeze-up  is  developed.  The  method  is  used  to  develop  a  computer  model  for  predicting  freeze-up 
in  the  Eisenhower  Lock  area  of  the  upper  St.  Lawrence  River  near  Massena.  The  results  show  that 
long-range  forecasts  of  freeze-up  in  a  river  can  be  made  with  good  accuracy. 

Problem  formulation 

For  a  well-mixed  river  the  conservation  of  thermal  energy  can  be  represented  by  a  one-dimen¬ 
sional  convection-diffusion  equation  (Brocard  and  Harleman  1976) 


!<»cp  o) 

where  t  »  time 

p  =  density  of  water 
Cp  *  specific  heat  of  water 
A  =  cross-sectional  area  of  the  river 
rw  *  water  temperature 
x  *  distance  along  the  river 
Q  =  river  discharge 

Ex  »  longitudinal  dispersion  coefficient 
B  »  channel  width 

<t>  *  net  heat  flux  from  the  river  per  unit  surface  area. 

If  we  assume  that  changes  in  river  discharge  are  not  significant,  we  can  simplify  eq  1  to  the  follow¬ 
ing  form: 


3  (E  — 
dx*  Xax  '  pCpD 


(2) 


where  U  is  the  average  flow  velocity  and  D  is  the  depth  of  flow.  By  neglecting  the  longitudinal 
mixing  term,  we  can  express  eq  2  in  its  Lagrangian  form: 


Dt  ~pCpD  • 


(3) 


For  a  lake  with  negligible  flow  velocity,  the  convective  term  can  be  neglected. 

The  surface  heat  exchange  is  a  complicated  function  of  ambient  atmospheric  conditions  (Paily 
et  al.  1974).  For  the  present  study  the  following  simple  relationship  is  assumed: 

♦  •Aw.O'w-r.)  W 


2 


W' 


Temperature  <*C) 


where  hv/k  is  an  energy  exchange  coefficient  and  7,  is  the  air  temperature.  Equations  3  and  4  de¬ 
fine  a  functional  relationship  between  the  water  temperature  and  the  ambient  air  temperature. 

For  long-range  freeze-up  forecasts,  an  analytical  description  of  the  predicted  air  temperature  is 
needed  for  eq  4.  Since  the  long-range  air  temperature  forecast  provided  by  the  National  Weather 
Service*  is  given  in  terms  of  deviations  from  the  normal  air  temperature,  an  analytical  representa¬ 
tion  of  the  normal  air  temperature  is  needed.  The  air  temperature  has  cyclic  variations  of  periods 
of  one  year  as  well  as  variations  with  shorter  periods.  This  cyclic  character  of  air  temperature  has 
been  well  recognized  (Kothandaraman  1971 ,  Song  et  al.  1973).  Time-dependent  variations  of  the 
air  temperature  can,  therefore,  be  represented  analytically  by  a  Fourier  series: 


_  N  2nirt  .  „  .  2nirt. 

r* =  T+„h  (Cn  C0ST 


(5) 


where  T  =  period  of  one  year 

T  =  mean  air  temperature 
Cn,Sn  =  Fourier  coefficients 

N  =  total  number  of  harmonics  in  the  series. 


In  this  study  the  time  variable  t  starts  on  1  October  each  year.  As  an  example,  Fourier  series 
approximations  with  N  =  1 0  for  both  the  air  temperature  and  the  water  temperature  records 
at  Massena,  New  York,  are  presented  in  Figure  1  for  the  year  between  1  October  1977  and  30 
September  1978.  Coefficients  of  the  Fourier  series  can  be  determined  from  the  temperature 
record  using  a  linear  regression  procedure  such  as  the  GLM  procedure  in  the  SAS  package  (SAS 
Institute  Inc.  1982). 


Figure  1.  Air  temperature  and  water 
temperature  at  Massena,  New  York, 

1  October  1977-30  September  1978. 


To  determine  the  appropriate  number  of  terms  to  be  included  in  the  Fourier  series  representa¬ 
tion  of  the  normal  air  temperature,  a  simple  statistical  analysis  car  be  made.  The  basis  of  the  air 
temperature  record  at  Massena,  New  York,  for  the  33  years  between  1948  and  1981 ,  a  Fourier 
series  of  N  =  12  is  obtained.  Table  1  shows  the  estimated  Fourier  coefficients  and  their  standard 
errors.  The  standard  errors  of  the  estimates  of  the  coefficients  of  harmonics  of  orders  higher  than 
one  are  of  the  same  order  as  the  estimated  values  of  the  coefficients  themselves.  This  indicates 
that  harmonics  of  periods  of  less  than  one  year  can  be  considered  to  be  random  components  and 
should  not  be  included  in  the  analytical  representation  of  the  normal  air  temperature.  To  further 
substantiate  this,  the  air  temperature  record  was  fitted  by  a  Fourier  series  of  orders  one  to  twelve. 

* Monthly  and  Seasonal  Weather  Outlook,  Climatic  Analyils  Center,  NOAA,  Washington,  DC. 
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Table  1.  Fourier  coefficients  for  air  tempera¬ 
ture  at  Massena,  hfew  York,  1948-1981. 

Standard  error 


Parameter  Estimate  _ of  estimate 


T 

6.63828048 

0.04126520 

Ci 

5.86978266 

0.05833369 

SI 

-13.09415854 

0.05838191 

C2 

0.16598747 

0.05833188 

S2 

0.58291816 

0.05838374 

C3 

-0.13276961 

0.05833649 

S3 

0.05808216 

0.05837906 

C4 

-0.29945449 

0.05835414 

S4 

0.29958279 

0.05836144 

CS 

-0.20162318 

0.05832959 

55 

-0.22805619 

0.05838598 

C6 

-0.00503364 

0.05833432 

56 

-0.09626112 

0.05838121 

a 

0.01613923 

0.05834101 

57 

-0.06291826 

0.05837450 

C8 

-0.17351080 

0.05836790 

58 

0.09971897 

0.05834760 

09 

0.05038423 

0.05834174 

59 

0.09057839 

0.05837374 

CIO 

-0.07903248 

0.05835328 

SlO 

0.22531156 

0.05836220 

Oil 

-0.00153682 

0.05836369 

511 

0.07895130 

0.05835277 

012 

0.13510373 

0.05936695 

512 

0.01369789 

0.05834843 

Table  2.  Accuracy  of  Fourier  series  repre¬ 
sentations  for  air  temperature  at  Massena, 
New  York,  1948-1981. 

Total  number  of  Correlation 

harmonics,  N _ parameter,  R* 


1 

0.8314 

2 

0.8329 

3 

0.8330 

4 

0.8337 

5 

0.8341 

6 

0.8341 

7 

0.8341 

8 

0.8343 

9 

0.8343 

10 

0.8345 

11 

0.8346 

12 

0.8346 

A  value  of  the  correlation  parameter  R 2  for  each  Fourier  series  was  calculated  (Table  2).  The 
parameter  R 2 ,  which  measures  how  much  variation  in  the  air  temperature  can  be  accounted  for 
by  a  Fourier  series,  is  the  ratio  of  the  sum  of  the  squares  of  the  predicted  value  divided  by  the 
sum  of  the  squares  of  the  recorded  values.  Table  2  shows  that  the  first  harmonic  accounts  for 
83%  of  the  variation  in  air  temperature;  the  inclusion  of  additional  harmonics  in  the  series  does 
not  improve  the  accuracy  of  the  representation  of  air  temperature  significantly.  Therefore,  the 
normal  air  temperature  can  best  be  represented  by  a  harmonic  function  with  a  period  of  one 
year.  Higher-order  variations  cannot  be  included  as  parts  of  the  normal  air  temperature.  Since 
the  higher-order  variations  change  from  year  to  year,  they  can  only  be  accounted  for  by  using 
forecasts  provided  by  the  National  Weather  Service. 

Analytical  treatment 

Using  only  the  first  harmonic  in  eq  5,  the  normal  air  temperature  variation  can  be  approximat¬ 
ed  by 


Th  =  T  +asin(~+0) 


(6) 


where  rN  =  periodic  function  representing  die  normal  temperature 

a  =  (Ci  +  Si  )** ,  the  amplitude  of  the  annual  temperature  cycle 
Q  =  tan'1  (Ci! Sx),  a  phase  angle. 

For  a  river  reach  with  known  water  temperature  at  the  upstream  end  (Fig.  2),  an  analytical  ex¬ 
pression  for  the  water  temperature  at  the  downstream  end  can  be  obtained  from  eq  3, 4  and  6; 

rj+i  -T\-a\  cosaj sin  +0i~  aj)  -*JfT 

- - - -  =  e  1  (7) 

Tj  o  -Tj-aj  cosaj  sin  (-j-  *  ) 

4 


. 


where  /  =  subscript  representing  the  location  of  the  river  reach 

Xj  and  Xj+ ,  =  space  coordinates  of  upstream  and  downstream  ends  of  the  reach 
fx.  =  travel  time  of  a  water  mass  In  the  reach,(xj+ ,  -  xj)/Ui 
Uj  =  mean  flow  velocity 
7j+i  =  water  temperature  at  (xj+  j ,  t) 

Ti  o  =  water  temperature  at  (xj,  t  -  rTj) 
kj  =  Awt/(pCpZ)j),  a  parameter  that  measures  the  sensitivity  of  water 
temperature  in  response  to  air  temperature  changes 
Z)j  =  mean  flow  depth 
Oj  =  tan-1  (In/k^T). 
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Figure  2.  Definition  sketch  of  a  river  reach. 

If  the  water  temperature  at  a  river  cross  section  is  known,  the  water  temperature  at  any  number 
of  downstream  stations  can  easily  be  obtained  by  successively  applying  eq  7.  For  a  river  system 
with  multiple  branches,  an  equation  for  conservation  of  thermal  energy  is  required  at  a  junction: 

To  =J,(7'“-*G“.*)/Gd  (8) 

water  temperature  at  the  outflow  side  of  the  junction 
total  number  of  incoming  brandies 

water  temperature  and  discharge  of  an  incoming  branch,  respectively 
outflow  from  the  junction. 

When  the  uppermost  station  is  at  the  outlet  of  a  lake,  the  following  solution  for  a  well-mixed 
lake  with  negligible  convective  velocity  applies: 

—  2irrL 

Tt  -  T  ~aL  cos aL  sin  ( ~jr~ +  0L-aL)  _  *-kL(tL- 

—  2irt0 

7*0 -r-<rL  cosaL  sin(—  +»L-aL) 

where  Th  *  lake  water  temperature  at  time  tL 
tQ  =  initial  time 

T0  *  lake  water  temperature  at  time  f0 

h-h„KfiCpDL) 

Dl  »  mean  depth  of  the  lake 


where 


= 
L  = 

Tujt  ®u,8  = 
QD  =  fiUl8  = 


5 


aL  =  tan"1  (2?r/kL  7) 

aL  and  0L  =  amplitude  and  phase  angle  of  the  air  temperature,  respectively. 

If  tQ  « tL ,  then  eq  9  becomes 

-  2irf. 

TL  =  T+aL  cos aL  sin (y  +  flL  -  aL).  (10) 

Equation  10  shows  that  the  amplitude  of  the  water  temperature  in  a  lake  is  equal  to  the  amplitude 
of  the  air  temperature  multiplied  by  cos  aL.  The  water  temperature  variation  lags  behind  the  air 
temperature  by  a  phase  angle  aL.  The  mean  water  temperature  is  equal  to  the  mean  air  tempera¬ 
ture.  When  the  headwater  of  a  river  is  not  a  lake,  one  can  use  eq  9  to  represent  the  water  tempera¬ 
ture  variation  at  the  headwater  by  treating  the  entire  drainage  basin  as  if  it  were  a  lake.  The  re¬ 
sponse  coefficient  kL  will  then  reflect  the  response  character  of  the  basin.  Song  et  al.  (1973) 
suggested  that  this  coefficient  decreases  as  the  basin  drainage  area  increases. 

Equations  7  and  9  describe  variations  of  water  temperature  when  the  air  temperature  is  expressed 
as  a  simple  harmonic  function  TNin  the  form  of  eq  6.  To  incorporate  variations  of  air  temperature 
with  a  period  less  than  1  year,  estimated  from  the  weather  forecast,  the  air  temperature  should  be 
approximated  by  a  combination  of  TN  and  a  series  of  short-term  deviations  as  shown  in  Figure  3. 


Figure  3.  Analytical  representation  of  a  series 
of  short-term  forecasts  for  air  temperature. 


The  short-term  deviation  from  TN  can  be  written  as 


6j ;  when  r,  <  t  <  tHl 

6Tl  = 

0  ;  otherwise . 

For  a  well-mixed  lake  the  effect  of  6  r„'  on  water  temperature  is 
f  0  ;  when  t  L  <  t{ 

5TL  =J6lll-e~*L"L~,,,1  ;when  </L  < 


01) 


(12) 


Sl[e"*l/,L“',+i U*Lf'L-V  j;  when  tL 


For  a  river  reach  as  defined  in  Figure  2,  the  effect  of  6  T*  is 


a)  when  t  -  fTj  <tt: 


fo  ;'<'i 

sai-e**^]  ;r ,<t<t 


i+i 


;t>t 


i+i 


b)  whenrt<f-fT 

fo 


;t<t-tT. 


87L.  =i8i[l-e'*j  <TJ  ]  ;  f  -  rT.  <  /  <  f4 


i+1 


;r>ri+1 


03) 


04) 


c)  when  1  “  fTj  >  *i+i : 

5rj+1  =  0.  (15) 

A  complete  solution  for  river  water  temperature  including  short-term  variations  can  be  obtained 
by  linear  superpositions  of  8Tj+1  to  Tj+1  and  6TL  to  T^. 


APPLICATION  TO  THE  UPPER  ST.  LAWRENCE  RIVER 

The  upper  St.  Lawrence  River,  which  is  the  outlet  from  Lake  Ontario,  extends  approximately 
160  km  from  Kingston,  Ontario,  to  the  Moses-Saunders  Power  Dam  near  Massena,  New  York 
(Fig.  4).  The  river  is  utilized  for  navigation  and  hydropower  production.  When  the  river  freezes 
up,  navigation  must  cease  and  head  losses  increase.  The  formation  of  the  ice  cover  also  reduces 
the  conveyance  capacity  of  the  river  and  affects  the  Lake  Ontario  water  level.  The  ability  to  pre¬ 
dict  the  freeze-up  date  is  therefore  important  in  managing  various  activities  on  the  river. 


Figure  4.  Plan  of  St.  Lawrence  River,  Lake  Ontario  to  Moses-Saunders  Power  Dam  (after  Wither¬ 
spoon  and  Poulin  1970). 
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Figure  5.  Schematic  representation  of  the  upper  St.  Lawrence  River. 

If  we  simplify  the  river  as  shown  in  Figure  5  and  assume  that  the  air  temperature  over  the 
entire  river  length  is  the  same,  then  the  normal  water  temperature  at  Masse na  can  be  obtained  from 
eq7: 

fR  N  =  r  +  a  cosaR  sin  (^?+  6  -  aR) 

+  exp (-*R/T)  [T’l.n  -  T-a  cos aR  sin  +  6- ctR) ]  (16) 

where  TR  N  =  normal  water  temperature  at  time  t 

tT  =  travel  time  between  Kingston  and  Massena,  (xm  -  x0)/U 
U  -  mean  flow  velocity 
*r  =KJ(pCpd) 

D  =  average  depth  of  flow  between  Kingston  and  Massena 
aR  =  tan'1  (2n/kRT) 


and  where 


rL  N  ~T  +  a  cos aL  sin (~  +  0  -aL) 


•  s2itta 


+  exp  [-*L(fL  " 'o)]  -T-a  cos aL  sin  (-y2  +  d  - aL)J 


07) 


where  T L  N  =  Kingston  normal  water  temperature  at  time  tL 

tL=t~tT 

ta  =  initial  time 

oL  =  tan'1  (2n/kLT) 

Dl  -  effective  depth  of  water  at  Kingston. 

Equations  16  and  17  are  analytical  expressions  for  water  temperature  at  Kingston  and  Massena 
when  the  air  temperature  is  approximated  by  To  include  the  effect  of  forecasted  deviations 

of  air  temperature  from  JN ,  the  Kingston  water  temperature  during  each  time  interval  can  be 
computed  by 

.  7’L-(T  +  8l)-«cosaLsin^'+d-aL) 

-*L  f'L-  _ _ L - 

...  _  2*f,  ' 

-(T*b^-acotah  «in(— ~-  +  6  -a^) 


(18) 
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where  Tt  =  Kingston  water  temperature  at  time  tL 
Tl*>  -  Kingston  water  temperature  at  time  r, 

6,  =  air  temperature  deviation  for  the  ith  time  period. 

The  water  temperature  at  Massena  can  be  approximated  by 

•  * 7 ■  — - Jr -  <”> 

Tl  -(T+  5,) -a  cos aR  sin  (~j~  +  6  - aR) 
when  <  /  <  fi+1 

where  7*R  is  the  Massena  water  temperature  at  time  t.  When  tt  is  in  the  (t'-l  )th  time  interval  and 
/  is  in  the  ith  time  interval,  eq  19  is  not  strictly  valid,  since  it  does  not  account  for  the  change  in 
air  temperature  deviation  from  one  time  interval  to  another.  However,  eq  19  is  used,  since  it  is 
easier  to  apply  and  the  accuracy  of  air  temperature  forecast  does  not  warrant  the  effort  of  calcu¬ 
lating  the  exact  solution  from  eq  13-15. 

Parameters  and  *R  can  either  be  estimated  from  the  average  depth  of  the  flow  in  the  river 
and  the  surface  heat  exchange  coefficient  or  be  determined  from  recorded  air  and  water  tempera¬ 
ture  data  by  nonlinear  regression  analyses.  The  latter  approach  is  more  accurate  and  is  used  in 
this  study.  Since  the  flow  depth  in  the  St.  Lawrence  River  remains  relatively  constant  over  time, 
kL  and  kR  values  are  considered  to  be  constants.  For  rivers  with  large  variations  of  flow  depth, 
these  coefficients  should  be  treated  as  functions  of  D.  Based  on  data  for  the  period  between  1965 
and  1981 ,  kL  and  kR  are  determined  by  the  NLIN  procedure  of  SAS  Institute  Inc.  (1982)  using 
half-month  intervals.  Table  3  summarizes  the  model  parameters  for  the  upper  St.  Lawrence  River. 


Table  3.  Model  parameters  for  the  upper  St.  Lawrence  River. 


Model  parameter 

Estimated  value 

Standard  error 

Mean  air  temperature  T 

6.639°C 

0.042°C 

Fourier  coefficient  C, 

5.869°C 

0.0S9°C 

Fourier  coefficient  S, 

-13.093°C 

0.059°C 

Response  parameter 

0.0191  day'1 

0.0010S  day'1 

Response  parameter  *R 

0.0345  day'1 

0.00182  day'1 

Previous  studies  (Joint  Board  of  Engineers  1926,  Witherspoon  and  Poulin  1 970,  Adams  1 976) 
have  shown  that  the  heat  exchange  coefficient  hwg  for  the  study  reach  is  about  21 .7  W/mIOC. 
Values  of  kL  and  kR  given  in  Table  3  correspond  to  an  effective  water  depth  of  23.56  m  at 
Kingston  and  an  average  depth  of  13.04  m  for  the  river  flow  between  Kingston  and  Massena. 

A  computer  program  for  forecasting  freeze-up  at  Massena  was  developed  and  is  presented  in 
Appendix  A.  The  freeze-up  at  Massena  is  defined  as  the  date  on  which  the  water  temperature  at 
the  Moses-Saunders  Dam  declined  to  0.3°  C  (Assel  1976, 1977).  At  that  time  it  is  probable  that 
an  ice  cover  has  formed  from  shore  to  shore  in  the  Eisenhower  Lock  area.  The  input  data  re¬ 
quired  for  the  freeze-up  forecast  include  the  initial  water  temperature  at  Kingston,  the  long-range 
forecast  for  air  temperature  deviation  from  the  normal  air  temperature  at  Massena,  and  the  pre¬ 
dicted  discharge  at  Moses-Saunders  Power  Dam.  The  air  temperature  deviation  can  be  estimated 
from  the  weather  outlook  provided  by  the  National  Weather  Service  along  with  the  standard  de¬ 
viation  of  air  temperature  given  in  Table  4.  Table  4  shows  large  fluctuations  in  Tt  during  December 
and  January. 


A  sample  freeze-up  forecast  beginning 
on  1  December  1983  is  included  in  Appen¬ 
dix  A.  The  measured  water  temperature 
at  Kingston  on  30  November  1983  was  TC, 
the  air  temperature  forecast  was  1 ,5°C  be¬ 
low  normal  for  1  December  to  IS  December 
1983  and  5°C  below  normal  for  16  Decem¬ 
ber  to  31  December,  and  the  forecasted 
travel  time  was  9  days.  The  model  fore¬ 
casted  a  freeze-up  date  of  24  December 
1983,  which  turned  out  to  be  the  day  of  the 
actual  freeze-up. 

To  test  the  accuracy  of  the  forecast  model, 
freeze-up  dates  were  calculated  for  the  16 
years  between  196S  and  1981  using  record¬ 
ed  air  and  water  temperature  data.  Simu¬ 
lated  and  recorded  water  temperature  de¬ 
clines  for  two  seasons  are  presented  in 
Figure  6.  The  observed  and  simulated 
freeze-up  dates  for  all  1 6  years  are  summar¬ 
ized  in  Table  S.  Confidence  limits  for  fore¬ 
casts  with  different  beginning  dates  of  simu¬ 
lation  are  presented  in  Table  6  and  Figure 
7.  Large  discrepancies  between  the  simu¬ 
lated  and  recorded  freeze-up  dates  during  a 


Table  4.  Semi-monthly  standard  deviations 
of  Maasena  air  temperature  from  first  har¬ 
monic. 


Stenderd 

Time  period _ Deviation  C  C) 


Oct. 

1-15 

4.S00 

16-31 

4.496 

Nov. 

1-1  s 

4.527 

16-30 

S.108 

Dec. 

1-1 5 

5.718 

16-31 

5.694 

Jan. 

1-1  s 

5.416 

16-31 

5.159 

Feb. 

1-15 

5.087 

16-29 

5.398 

Mar. 

1-1 S 

4.973 

16-31 

4.854 

Apr. 

1-1 S 

4.097 

16-30 

4.121 

May 

1-1 S 

4.471 

16-31 

4.361 

Jun. 

1-1S 

3.900 

16-30 

3.440 

Jul. 

1-15 

3.195 

16-31 

3.085 

Aug. 

1-15 

3.307 

16-31 

3.691 

Sep. 

1-15 

4.021 

16-30 

4.246 

a.  1966-67. 


b.  1978-79. 


Figure  6.  Recorded  and  simulated  water  temperature,  Massena,  New  York. 
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eMmp  dates,  St.  Lawrence  River  at  Massetu,  1 965-1980. 


Observed 
frees e-up 

Forecasted  freere-up 

Tur 

1  Oct 

16  Oct 

1  Nov 

16  Nov 

1  Dec 

16  Dec 

1965-66 

8  Jen 

1  Jen 
-7 

1  Jen 
-7 

1  Jen 
-7 

1  Jen 
-7 

1  Jen 
-7 

3  Jan 
-5 

1966-67 

30  Dee 

29  Dec 
-1 

29  Dec 
-1 

30  Dec 

0 

29  Dec 
-1 

30  Jan 

0 

11  Jen 
+12 

1967-68 

3  Jan 

29  Dec 
-5 

28  Dec 
-6 

28  Dec 
-6 

27  Dec 
-7 

29  Jan 
-5 

1  Jan 
-2 

1968-69 

25  Dec 

28  Dec 
+3 

26  Dec 
+1 

26  Dec 
+1 

25  Dec 

0 

26  Dec 
+1 

26  Dec 
+1 

1969-70 

28  Dec 

25  Dec 
-3 

24  Dec 
-4 

23  Dec 
-5 

23  Dec 
-5 

23  Dec 
-5 

27  Dec 
-5 

1970-71 

23  Dec 

30  Dec 
+7 

30  Dec 
+7 

29  Dec 
+6 

29  Dec 
+6 

25  Dec 
+2 

27  Dec 
+4 

1971-72 

8  Jan 

31  Jen 
-8 

30  Dec 
-9 

31  Jen 
-8 

30  Dec 
-9 

31  Jan 
-8 

8  Jan 

0 

1972-73 

21  Dec 

24  Dec 
+3 

23  Dec 
+2 

23  Dec 
+2 

22  Dec 
+1 

24  Dec 
+3 

25  Dec 
+4 

1973-74 

30  Dec 

1  Jen 
+2 

1  Jen 
+2 

1  Jen 
+2 

30  Dec 
0 

1  Jan 
+2 

1  Jan 
+2 

1974-75 

22  Jen 

18  Jen 
-4 

14  Jan 
-8 

14  Jen 
-8 

10  Jen 
-12 

21  Jan 
-1 

23  Jan 
+1 

1975-76 

27  Dec 

27  Dec 
0 

27  Dec 
0 

27  Dec 

0 

28  Dec 
+1 

27  Dec 
0 

28  Dec 
+1 

1976-77 

14  Dec 

22  Dec 
+8 

21  Dec 
+7 

17  Dec 
+3 

18  Dec 
+4 

15  Dec 
+1 

15  Dec 
+1 

1977-78 

30  Dec 

1  Jen 
+2 

29  Dec 
+1 

30  Dec 

0 

28  Dec 
-2 

No  Kingston 
Data  11/23-12/31 

1978-79 

28  Dec 

28  Dec 
-2 

27  Dec 
-3 

27  Dec 
-3 

27  Dec 
-3 

26  Dec 
-4 

10  Jan 
+1 

1979-80 

12  Jen 

12  Jen 
0 

11  Jan 
-1 

11  Jan 
-1 

12  Jen 
0 

13  Jan 
+1 

10  Jan 
-2 

1980-81 

21  Dec 

16  Dec 
-5 

16  Dec 
-5 

16  Dec 
-5 

16  Dec 
-5 

16  Dec 
-5 

16  Dec 
-5 

Mean 

-0.625 

-0.625 

-1.813 

-2 .438 

-1 .667 

+0.800 
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Table  6.  Confidence  limits  of  the  freeze-up  forecast. 


Oct  1 

Oct.  16 

Nov.  1 

Nov.  16 

Dec.  1 

Dec.  16 

95%  upper  limit 

1.86 

0.94 

0.47 

0.12 

0.34 

3.07 

80%  upper  limit 

0.94 

-0.01 

-0.38 

-0.83 

-0.41 

2.22 

50%  upper  limit 

0.18 

-0.79 

-1.07 

-1.61 

-1.02 

1.53 

Mean  error 

-0.63 

-1.63 

-1.81 

-2.44 

-1.67 

0.80 

50%  lower  limit 

-1.43 

-2.46 

-2.55 

-3.27 

-2.31 

0.07 

80%  lower  limit 

-2.19 

-3.24 

-3.25 

-4.05 

-2.92 

-0.62 

95%  lower  limit 

-3.11 

-4.19 

-4.09 

-5.00 

-3.67 

-1.47 

Figure  7.  Confidence  limits  for  the 
deviation  of  forecasted  freeze-up 
date  from  the  observed  date. 


few  winters  may  be  attributed  to  the  simplicity  of  the  freeze-up  criteria  used.  The  water  tempera¬ 
ture  at  Massena  may  fluctuate  around  0.3°C  over  a  period  of  several  weeks.  During  this  period  the 
forecast  area  may  alternately  freeze  and  thaw.  The  use  of  a  single  date  for  the  freeze-up  is  inade¬ 
quate  in  describing  this  phenomenon.  The  inability  of  a  constant  surface  heat  exchange  coefficient 
to  account  for  some  of  the  meteorological  phenomena,  such  as  snowfalls  and  extreme  wind  condi¬ 
tions,  may  also  contribute  to  the  discrepancies  in  the  freeze-up  forecast. 


SUMMARY 

In  this  report  a  method  for  making  long-range  forecasts  of  river  freeze-up  is  developed.  It  is 
shown  that  the  water  temperature  decline  in  a  river  is  governed  by  the  convection  of  the  water 
mass  in  the  river  and  the  heat  exchange  at  the  free  surface.  The  surface  heat  exchange  can  be  ex¬ 
pressed  in  terms  of  the  difference  between  the  water  temperature  and  the  ambient  air  temperature. 
The  air  temperature  can  be  effectively  represented  by  a  linear  combination  of  the  normal  air  tem¬ 
perature  described  by  a  harmonic  function  with  a  period  of  one  year  and  short-term  variations  ob¬ 
tained  from  the  National  Weather  Service.  The  application  of  this  mehtod  to  the  St.  Lawrence  River 
between  Kingston  and  Massena  shows  that  it  is  capable  of  providing  accurate  freeze-up  forecasts. 
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APPENDIX  A.  BASIC  PROGRAM  FOR  ST.  LAWRENCE  RIVER  FREEZE-UP  FORECAST 


LIST 

1  REM  M******t*************************»*«*******************««**«* 

2  REM 

3  REM  INPUT  TO  ST  LAWRENCE  RIVER  FREEZE-UP  MODEL 

4  REM  TO  -  INITIAL  DATE  OF  FORECAST  IN  DAYS  FROM  OCTOBER  1 

5  REM  D1-D8  -  FORECASTED  AIR  TEMP  DEVIATION  FROM  NORMAL  IN  DEG  C 

6  REM  TWO  -  KINGSTON  WATER  TEMP  IN  DEG  C  ON  FIRST  DAY  OF  FORECAST 

7  REM  TT  -  TRAVEL  TIME  IN  DAYS  FROM  KINGSTON  TO  MASSENA  <2. 1303E6/Q) 

8  REM  IYEAR7.  -  YEAR  OF  FORECAST 

9  REM  ************************************************************** 

10  KL  -  0.01910924 

15  TIL  «  2.72018  -  ATN  (0.01721  /  KL> 

20  KR  -  0.03452020 

25  T2R  -  2.72018  -  ATN  (0.01721  /  KR) 

30  CL  -  1  /  SQR  (1  +  <0.01721  /  KL>  ~  2) 

40  CR  -  1  /  SQR  <1  <0.01721  /  KR)  ~  2) 

44  REM 

45  REM  INPUT  TEMP  DEVIATIONS  AND  INITIAL  CONDITIONS 

46  REM 

50  READ  T0,D1,D2,D3,D4,D5,D6,D7,D8,TW0,TT, I YEAR* 

60  IF  TO  >  -  1  AND  TO  <  -  31  THEN  ID%  -  TO 

70  IF  TO  >  «  32  AND  TO  <  -  61  THEN  ID7.  »  TO  -  31 

80  IF  TO  >  -  62  AND  TO  <  -  92  THEN  ID 7.  -  TO  -  61 

90  IF  TO  >  -  93  THEN  1DX  -  TO  -  92 

94  PRINT 

95  PRINT 

96  PRINT 

100  IF  TO  >  -  1  AND  TO  <  -  31  THEN  PRINT  "THE  OCTOBER  "jID7.j 

110  IF  TO  >  -  32  AND  TO  <  -  61  THEN  PRINT  "THE  NOVEMBER  “j  ID7.| 

120  IF  TO  >  -  62  AND  TO  <  -  92  THEN  PRINT  "THE  DECEMBER  "jID7.j 

130  IF  TO  >  -  93  THEN  PRINT  "THE  JANUARY  ID7.J 

140  PRINT  "  FORECAST  AT  MASSENA,  N.Y.  FOR  "jIYEARY. 

150  IB  -  INT  (TO) 

154  REM 

155  REM  START  LOOP  TO  COMPUTE  MASSENA  WATER  TEMPERATURE 

156  REM 

160  FOR  IT  -  IB  TO  123 

170  TIME  ■  IT 

180  TMTT  -  TIME  -  TT 

184  REM 

185  REM  UPDATE  LAKE  TEMP,  AIR  TEMP  DEVIATIONS,  AND  CHECK  FOR  NEW  TO 

186  REM 


190 

IF 

TIME 

-  16 

OR 

TIME  -  32 

OR 

TIME  -  47  THEN  GOSUB 

660 

200 

IF 

TIME 

-  62 

OR 

TIME  -  77 

OR 

TIME  -  93  THEN  GOSUB 

660 

210 

IF 

TIME 

>  - 

1  AND  TIME  < 

m 

15  THEN  DR  -  D1 

220 

IF 

TIME 

>  - 

16 

AND  TIME  < 

m 

31  THEN  DR  -  D2 

230 

IF 

TIME 

>  - 

32 

AND  TIME  < 

m 

46  THEN  DR  -  D3 

240 

IF 

TIME 

>  - 

47 

AND  TIME  < 

m 

61  THEN  DR  »  D4 

250 

IF 

TIME 

>  - 

62 

AND  TIME  < 

m 

76  THEN  DR  -  D5 

260 

IF 

TIME 

>  - 

77 

AND  TIME  < 

* 

92  THEN  DR  -  D6 

270 

IF 

TIME 

>  - 

93 

AND  TIME  < 

m 

107  THEN  DR  -  D7 

280 

IF 

TIME 

>  107  THEN  DR  *  D8 

290 

IF 

TMTT 

<  - 

15 

THEN  DL  - 

D1 

300 

IF 

TMTT 

>  - 

16 

AND  TMTT  < 

m 

31  THEN  DL  -  D2 

310 

IF 

TMTT 

>  - 

32 

AND  TMTT  < 

m 

46  THEN  DL  -  D3 

320 

IF 

TMTT 

>  - 

47 

AND  TMTT  < 

a> 

61  THEN  DL  -  D4 

330 

IF 

TMTT 

>  • 

62 

AND  TMTT  < 

* 

76  THEN  DL  -  D5 

IS 


340 

IF 

TMTT 

> 

-  77  AND  TMTT  < 

-  92  THEN  DL  - 

06 

350 

IF 

TMTT 

> 

-  93  AND  TMTT  < 

-  107  THEN  DL 

-  07 

360 

IF 

TMTT 

> 

107  THEN  DL  -  08 

370 

IF 

TIME 

81 

16  THEN  TO  -  15 

380 

IF 

TIME 

m 

32  THEN  TO  -  31 

390 

IF 

TIME 

m 

47  THEN  TO  -  46 

400 

IF 

TIME 

m 

62  THEN  TO  -  61 

410 

IF 

TIME 

m 

77  THEN  TO  -  76 

420 

IF 

TIME 

m 

93  THEN  TO  -  92 

430 

IF 

TIME 

m 

107  THEN  TO  -  106 

434  REM 

433  REM  COMPUTE  MASSENA  RIVER  TEMP 
436  REM 

440  A1  -  0.01721  *  TIME  +  T2R 
430  A3  *  0.01721  *  TMTT  +  T2R 
460  ER  -  EXP  <  -  KR  *  TT) 

470  XD  »  TIME 

480  GOSUB  700 

490  A  *  6.639  DR 

300  B  -  <TL  -  6.639  -  DR)  *  ER 

310  C  -  14.348  *  CR  *  SIN  (Al) 

320  D  ■  14.348  *  ER  *  CR  *  SIN  (A3) 

330  TR-A  +  B  +  C-  D 

335  PRINT  "  TIME*  "jTIMEs ”  DAYS, RIVER  TEMP-  " ;  TR; "  DEG  C" 

340  IF  TR  <  *  0.3  THEN  GOTO  370 

330  NEXT  IT 

360  PRINT  -THE  DATE  OF  FREEZE-UP  WAS  NOT  REACHED  BY  JANUARY  31 
364  REM 

563  REM  DETERMINE  CALENDAR  DATE  OF  FREEZE-UP 
366  REM 

370  IF  TIME  <  -  92  GOTO  610 

380  ID7.  -  TIME  -  92 

390  PRINT  -THE  DATE  OF  FREEZE-UP  IS  JANUARY  "|ID*>t 

600  GOTO  630 

610  ID7.  -  TIME  -  61 

620  PRINT  “THE  DATE  OF  FREEZE-UP  IS  DECEMBER  ”;ID7. 

630  REM 
640  GOTO  SO 
630  END 
634  REM 

633  REM  SUBROUTINE  TO  UPDATE  LAKE  TEMP  AT  EACH 

636  REM  CHANGE  OF  AIR  TEMP  DEVIATION 

637  REM 

660  XD  -  TIME  ♦  TT 
670  GOSUB  700 
680  TWO  -  TL 
690  RETURN 
694  REM 

693  REM  SUBROUTINE  TO  COMPUTE  LAKE  TEMP 
696  REM 

700  A2  -  0.01721  *  <XD  -  TT)  ♦  TIL 

710  A4  •  0.01721  *  <T©>  ♦  TIL 

720  EL  -  EXP  <  -  KL  *  <XD  -  TT  -  TO)) 

730  A  -  6.639  «•  DL 

740  B  -  (TWO  -  6.639  -  DL)  *  EL 

750  C  ■  14.348  *  CL  *  SIN  <A2) 

760  D  -  14.34^  *  EL  *  CL  *  SIN  (A4) 

770  TL-A-*B*C-D 
780  RETURN 

790  DATA  62,0, 0,0,0,-1.3,-3,0,0,7,9.8384 
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RUN 


THE  DECEMBER  1  FORECAST  AT  MASSENA,  N.Y.  FOR  8384 


TIME* 

62 

DAYS, RIVER 

TEMP- 

5. 3980055 

DEQ  l 

C 

TIME¬ 

63 

DAYS, RIVER 

TEMP- 

5.22429169 

DEG 

c 

TIME* 

64 

DAYS, RIVER 

TEMP- 

5.0504605 

DEG  l 

2 

TIME* 

65 

DAYS, RIVER 

TEMP- 

4.8765714 

DEG  1 

C 

TIME¬ 

66 

DAYS, RIVER 

TEMP- 

4.70268361 

DEG 

c 

TIME* 

67 

DAYS, RIVER 

TEMP- 

4.52885622 

DEG 

c 

TIME* 

68 

DAYS, RIVER 

TENP- 

4.35514823 

DEG 

c 

TIME¬ 

69 

DAYS, RIVER 

TEMP- 

4. 18161835 

DEG 

c 

TIME¬ 

70 

DAYS,  RIVER- 

TEMP- 

4 . 00832524 

DEG 

c 

TIME* 

71 

DAYS.  RIVER 

TEMP- 

3.81451759 

DEG 

c 

TIME¬ 

72 

DAYS, RIVER 

TEMP- 

3. 62145706 

DEG 

c 

TIME* 

73 

DAYS.  RIVER- 

TEMP- 

3. 42919415 

DEG 

c 

TIME- 

74 

DAYS.  RIVER 

TEMP- 

3.23777917 

DEG 

c 

TIME¬ 

75 

DAYS,  RIVER- 

TEMP- 

3.0472624  ! 

DEG  ( 

2 

TIME¬ 

76 

DAYS,  RIVER 

TEMP- 

2.8576939  ] 

DEG  ( 

2 

TIME* 

77 

DAYS,  RIVER- 

TEMP- 

1 . 57073723 

DEG 

c 

TIME- 

78 

DAYS,  RIVER 

TEMP- 

1.38631337 

DEG 

c 

TIME- 

79 

DAYS, RIVER 

TEMP- 

1 . 20292832 

DEG 

c 

TIME- 

80 

DAYS, RIVER 

TEMP- 

1 . 02063256 

DEG 

c 

TIME- 

81 

DAYS, RIVER 

TEMP- 

.839476316 

DEG 

c 

TIME* 

82 

DAYS, RIVER 

TEMP- 

. 659509452 

DEG 

c 

TIME* 

83 

DAYS, RIVER 

TEMP- 

.480781701 

DEG 

c 

TIME- 

84 

DAYS,  RIVER- 

TEMP- 

. 303342351 

DEG 

c 

TIME- 

85 

DAYS,  RIVER 

TEMP- 

. 127240472 

DEG 

c 

THE  DATE  OF  FREEZE-UP  IS  DECEMBER  24 
?OUT  OF  DATA  ERROR  IN  50 
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A  facsimile  catalog  card  in  Library  of  Congress  MARC 
format  is  reproduced  below. 


Hung  Tao  Shen 

Forecasting  water  temperature  decline  and  freeze-up 
in  rivers  /  by  Hung  Tao  Shen,  E.P.  Foltyn  and  S.F. 

Daly.  Hanover,  N.H.:  Cold  Regions  Research  and  Engi¬ 
neering  Laboratory;  Springfield,  Va. :  available  from 
National  Technical  Information  Service,  1984, 
iii,  22  p. ,  illus.;  28  cm.  (CRREL  Report  84-19.) 
Bibliography:  p.  12. 

1.  Ice  formation.  2.  Rivers.  3.  River  ice. 

4.  St,  Lawrence  River.  I.  Foltyn,  E.P.  II.  Daly, 

5. F.  ITL  United  States.  Army.  Corps  of  Engineers. 

IV.  Cold  Regions  Research  and  Engineering  Laboratory, 
Hanover,  New  Hampshire.  V.  Series:  CRREL  Report  84-19. 
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